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Abstract. We study by Monte Carlo simulations the influence of bond dilution on the three-dimensional 
Ising model. This paradigmatic model in its pure version displays a second-order phase transition with 
a positive specific heat critical exponent a. According to the Harris criterion disorder should hence lead 
' O ' to a new fixed point characterized by new critical exponents. We have determined the phase diagram 

C ' of the diluted model, between the pure model limit and the percolation threshold. For the estimation 

of critical exponents, we have first performed a finite-size scaling study, where we concentrated on three 
different dilutions. We emphasize in this work the great infiuence of the cross-over phenomena between 
the pure, disorder and percolation fixed points which lead to effective critical exponents dependent on the 
concentration. In a second set of simulations, the temperature behaviour of physical quantities has been 
^ ' studied in order to characterize the disorder fixed point more accurately. In particular this allowed us to 

' estimate ratios of some critical amplitudes. In accord with previous observations for other models this 

provides stronger evidence for the existence of the disorder fixed point since the amplitude ratios are more 
sensitive to the universality class than the critical exponents. Moreover, the question of non-self-averaging 
at the disorder fixed point is investigated and compared with recent results for the bond-diluted q — 4 
Potts model. Overall our numerical results provide evidence that, as expected on theoretical grounds, the 
critical behaviour of the bond-diluted model is governed by the same universality class as the site-diluted 
model. 
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1 Introduction the two-dimensional (2D) Ising model has attracted great 

interest in the past years because of its intermediate situ- 
, The influence of quenched, random disorder on phase tran- ation (a = 0) which implies a marginal influence of disor- 
b : sitions has been the subject of numerous experimental and ^er g] . From an experimental pomt of view, a confirma- 
■ theoretical investigations since more than 20 years ffl. ^lon of the Harris criterion was reported in a LEED mves- 
They concern both first- and second-order phase transi- tigation of a 2D order-disorder transition In the case 
tions, especially in two dimensions. The qualitative influ- ^ first-order phase transition m the pure model, disorder 
ence of quenched, short-range correlated random dis- expected to soften the transition and, under some cn- 
order at second-order phase transitions is well understood cumstances, may even induce a second-order transition |7| . 
since Harris 3 proposed a relevance criterion based on For 2D systems, the latter scenario has been proved by 
the knowledge of the specific heat critical exponent a of Aizenman and Wehr on rigorous theoretical grounds 
the pure model: when a is positive, under a coarse grain- M test these theoretical predictions also the Potts 
ing, the disordered system should reach a so-called finite- "^^del HH has been intensively studied m 2D [U since 
randomness disorder flxed point characterized by altered displays the two different regimes: a second-order phase 
critical exponents, whereas if a is negative, the univer- transition when the number of states per spin g < 4 and 
sality class of the pure system will persist. In particular ^ first-order one when q > 4. These results were obtained 
by different techniques including Monte Carlo simulations, 
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theory, and high-temperature series expansions |12II18I 

In three dimensions (3D), the disordered Potts model 
has of course been studied only later: the case q = 3, corre- 
sponding to a very weak first-order transition for the pure 
model (23) . has been investigated numerically by Balles- 
teros et al. for site dilution, and the case 5 = 4, 
exhibiting in the pure system a strong first-order tran- 
sition has been studied very recently by us [2323 
I28lf? 9 , 30 for bond dilution via large-scale Monte Carlo 
simulations. Only the 3D Ising model with site dilution 
has also been extensively studied using Monte Carlo sim- 
ulations ^;n'^TTf^m^''s^wM'W'w:^'^2'w^'^^ , 

or field theoretic renormalization group approaches 4fii 

treated in the low-dilution regime (concentration of mag- 
netic bonds, p, close to 1) by analytical perturbative renor- 
malization group methods |56ll57ti 58 ,59 ,60 where a new 
fixed point independent of the dilution has been found, 
but for stronger disorder only Monte Carlo simulations 
remain valid. 

The first numerical studies p^33ll34| suggested a con- 
tinuous variation of the critical exponents along the crit- 
ical line but after the works 35,36,37' 38"," 39| it became 
clear that the concentration dependent critical exponents 
found in Monte Carlo simulations are effective ones, char- 
acterizing the approach to the asymptotic regime. The 
critical exponents (3 and 7 associated with the magnetisa- 
tion and susceptibility, respectively, were shown by Heuer |38| 
to be concentration dependent in the region 0.5 < p < 1- 
The conclusion of Heuer was: while a crossover between 
the pure and weakly random fixed points accounts for the 
behaviour of systems above p ~ 0.8, in more strongly dis- 
ordered systems a more refined analysis is needed. We 
should mention here that the meaning of p is not exactly 
the same in the papers mentioned above which refer to site 
dilution and in our present study in which we are inter- 
ested in bond dilution. That is why the numerical values 
of p given before for the different regimes cannot directly 
be taken over to the bond-dilution case. 

Another important question has been investigated by 
Wiseman and Domany |43ll44j : it concerns the question 
of the possible lack of self-averaging which can happen in 
disordered systems. For the 3D site-diluted Ising model 
close to criticality they explicitly showed that physical 
quantities such as the magnetisation or susceptibility are 
indeed not self-averaging. Although simulations ^| re- 
vealed that disorder realized in a canonical manner (fixing 
the fraction p of magnetic sites) leads to different results 
than those obtained from disorder realized in a grand- 
canonical ensemble (assigning to each site a magnetic mo- 
ment with probability p), the renormalization-group ap- 
proach of Ref. |HJ shows that the canonical constraint 
is irrelevant, even near the random fixed point, suggest- 
ing that the observed differences are a finite-size effect. 
The studies of Ref. 03 were based on the crucial observa- 
tion that it is important to take into account the leading 
corrections-to-scaling term in the infinite-volume extrap- 
olation of the Monte Carlo data. Thus, the main problem 



encountered in these studies of the disordered Ising model 
was the question of measuring effective or asymptotic ex- 
ponents. Although the change of universality class should 
happen theoretically for arbitrarily weak disorder, the new 
critical exponents appear only in a small temperature re- 
gion around the critical point, whose size is controlled by 
the concentration of the non-magnetic compound. Equiv- 
alently, in finite-size scaling studies very large lattices are 
required to observe the asymptotic behaviour. In fact, the 
asymptotic regions cannot always be reached practically 
and one therefore often measures only effective exponents. 

Another crucial problem of the new critical (effective 
or asymptotic) exponents obtained in these studies is that 
the ratios (3/v and j/v occuring generically in finite-size 
scaling analyses are almost identical for the disordered and 
pure models. For the pure 3D Ising model, accurate values 
are 

V = 0.6304(13), 7? = 0.0335(25), 

which gives: 

= 0.517(3), 7/j/ = 2-ry = 1.966(3), 

and OL = 0.1103(1) > 0, i.e., disorder should be relevant 
according to the Harris criterion. 

For the site-diluted model, the asymptotic exponents 
given by Ballesteros et al. 021 are: 

= 0.519(3), 7/1/ = 1.963(5), = 0.6837(53). 

Thus, finite-size scaling techniques will only be able 
to differentiate between the values of v at the two fixed 
points, but will not be very efficient for distinguishing ra- 
tios of critical exponents. Even if /3 and 7 themselves hap- 
pen to be quite different, the ratios fijv and ^ jv are very 
close. That is the reason why a study of the temperature 
behaviour of the magnetisation and susceptibility will be 
very helpful for an independent determination of the ex- 
ponents /3 and 7. 

Contrary to previous studies of the disordered Ising 
model which were concerned with site dilution, we have 
chosen to model the disorder by bond dilution in order 
to compare these two kinds of disorder and to verify that 
they indeed lead to the same set of new critical expo- 
nents, as expected theoretically by universality arguments. 
In the following we shall thus consider the bond-diluted 
Ising model in 3D whose Hamiltonian with uncorrelated 
quenched random interactions can be written (in a Potts 
model normalization) as 

(ij) 

where the spins take the values (7^ = ±1 and the sum 
goes over all nearest- neighbour pairs («, j). The coupling 
strengths are allowed to take (grand-canonically) two dif- 
ferent values Kij = K = J /ksT and with probabilities 
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p and I — p, respectively, 

nK^,] - n ^(^^^■) 

= n [^"^(^'J- + p)5{K,,)], (2) 

c = I —p being the concentration of missing bonds, which 
play the role of the non-magnetic impurities. 

The plan of the rest of the paper is as follows: in sec- 
tion 2, we present the phase diagram of the bond-diluted 
Ising model. Then, in section 3, we discuss the averaging 
procedure and investigate the question of possible non- 
self-averaging of physical quantities. The section 4 is de- 
voted to the critical behaviour of the disordered model. 
In the first part we present our finite-size scaling study 
of three particular concentrations p, and the second part 
deals with the temperature scaling for the same three di- 
lutions. Finally, section 5 contains our conclusions. 



2 Phase diagram 

In order to determine the phase diagram and the criti- 
cal properties at a few selected dilutions we performed in 
this study large-scale Monte Carlo simulations on sim- 
ple cubic lattices with V = spins (up to L = 96) 
and periodic boundary conditions in the three space di- 
rections, using the Swendsen-Wang cluster algorithm [HSI 
for updating the spins. The histogram reweighting tech- 
nique )64II65| was employed to extend the results over 
a range of K around the simulation point. All physical 
quantities are averaged over 2 000 - 5 000 disorder real- 
isations, indicated by a bar (e.g., fh for the magnetisa- 
tion). Standard definitions were used, e.g., for a given 
disorder realisation, the magnetisation is defined accord- 
ing to m = where (...) stands for the thermal av- 
erage and /i = {N^ - Ni)/{N'^ + Ni). The susceptibil- 
ity follows from the fluctuation-dissipation relation, x — 
KV{{fj?) — (|/i|)^). Here, one should notice that a defini- 
tion which is usually used in the disordered phase in Monte 
Carlo simulations, and which seems to be more stable, 
X = KV{^'^), should be avoided in quenched disordered 
systems. Simulating at a given temperature above the crit- 
ical point, one may indeed encounter samples which have 
higher effective transition temperatures, and which thus 
have non-zero (|/i|). 

The phase diagram was obtained by locating the max- 
ima of the average susceptibility xl (a diverging quantity 
in the thermodynamic limit) for increasing lattice sizes L 
as a function of the coupling strength K, with the dilution 
parameter p varying from the neighbourhood of the pure 
model {p — 0.95) to the very diluted model {p — 0.36), see 
Fig. ID 

Below the percolation threshold pc ~ 0.2488 jSSj, one 
does not expect any finite-temperature phase transition 
since without any percolating cluster in the system long- 
range order is impossible. The approximate phase dia- 
gram |67j as obtained from the susceptibility maxima for 



the largest lattice size {L ~ 20) is shown in Fig. [3 For 
comparison, we have drawn a simple mean-field (MF) es- 
timate of the transition point 



ifr(p)=p^c(i), 



(3) 



where Kc{l) = 0.443 308 8(6) [SHI is the accurately known 
transition point of the pure model (in the Potts model nor- 
malization), and the single-bond effective-medium (EM) 
approximation (69j 



(l-p,)e^^(i)-(l-p) 



P-Pc 



(4) 



which gives a very good agreement with the simulated 
transition line over the full dilution range. We have omit- 
ted results from recent high-temperature series expansions 
[70] since, on the scale of Fig. [3 they would just fall on 
top of the Monte Carlo data. 

To get an accurate determination of Kc{L), we used 
the histogram reweighting technique with at least A^mcs = 
2500 Monte Carlo sweeps and between 2000 and 5000 sam- 
ples of disorder. The number of Monte Carlo sweeps is jus- 
tified by the increasing behaviour of the energy autocor- 
relation time Tg as a function of p and L. For each size, we 
performed at least 250 independent measurements of the 
physical quantities (A^mcs > 250 Te). For a second-order 
phase transition, the autocorrelation time is expected to 
behave as at the critical point where z is the dynam- 
ical critical exponent. For the disordered Ising model, we 
obtained from the least-squares fits shown in Fig. O the 
values of z compiled in Tabled We see that the critical 
slowing-down weakens for the disordered model and that 
z becomes effectively smaller when the concentration of 
magnetic bonds p decreases. The observed variation of the 
dynamical exponent z is probably due to the influence of 
the different flxed points encountered (pure, disorder and 
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Fig. 1. Variation of the average magnetic susceptibility xl 
versus the coupling strength K — J/ksT for several concen- 
trations p = 0.95, 0.90, . . . , 0.36 and L = 10, 12, 14, 16, 18, 20. 
For each value of p and each lattice size, the curves are obtained 
by standard histogram reweighting of the simulation data at 
one value of K. 
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Fig. 2. Phase diagram of the 3D bond-diluted Ising model 
compared with the mean-field and effective-medium approxi- 
mations. The open circle marks the location of the percolation 
threshold. 



percolation), as will be discussed below in more detail for 
the static critical exponents. 

The largest autocorrelation time observed for the dis- 
ordered model was around ~ 9 for p = 0.7 and L = 96. 



10 




10 



100 



Fig. 3. Energy autocorrelation time re versus the size L of the 
system on a log- log scale. The pure case corresponds to p = 1. 
The straight lines show fits of the form Te oc L" , yielding the 
effective dynamical critical exponents z compiled in Tabled 



3 Non-self-averaging 

In order to achieve accurate results for quenched, disor- 
dered systems in numerical simulations it is important to 
obtain an estimate of the required number of disorder re- 
alisations. This is particularly important in the vicinity of 
a critical point where the correlation length diverges. As 
a consequence the (disorder) distributions of physical ob- 
servables typically do not become sharper with increasing 
system size at a finite-randomness disorder fixed point. 



Table 1. Effective dynamical critical exponent z as obtained 
from linear fits of log vs log L. The simulations are performed 
at the pseudo-critical couplings Kc{L) (cf. Table|21and Fig. llUI 
below). Here we give in the second line only the simulation 
coupling K for the largest lattice size L = 96. 



p 


1 


0.7 


0.55 


0.4 


K 


0.4433 


0.6535 


0.8649 


1.3136 


z 


0.59 


0.41 


0.38 


0.27 
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Fig. 4. Disorder distribution of the susceptibility for the Ising 
and q = A Potts models with a concentration of magnetic bonds 
of 0.55 and 0.56, respectively, and a lattice of size L — 64. The 
simulations are made ai K — 0.8655 « Kc{L) for the Ising 
model and K = 1.12945 « Kc{L) for the Potts model. In the 
latter case, the concentration p — 0.56 belongs to the second- 
order regime 126 j . The running average over the samples Xj is 
shown by the thick solid line. 



Rather their relative widths stay constant, a phenomenon 
called non-self-averaging. In order to investigate the disor- 
der averages, we produced Ns different samples and com- 
puted the corresponding susceptibilities Xj j 1 ^ J ^ ^s- 
In Fig. ^ we compare the distributions for the Ising and 
q = 4 Potts models (in the second-order regime |5H]). The 
figure shows that the dispersion of the values of x is less 
important for the disordered Ising model where rare events 
correspond to low values of x- This implies that their con- 
tribution to the average is not so crucial as for the Potts 
model where a long rare-event tail is found on the large-x 
side. This is the reason why for the Ising model the fluctu- 
ations in the average value disappear after a few hundred 
realisations and why the disorder averaging procedure is 
more efficient than for the Potts model. 

To test if self-averaging is valid or not for the dis- 
ordered model |43II44| . we investigated the variation of 
the shape of the probability distributions of the suscep- 
tibility when the size is increasing. Dividing the interval 
(Xmax — Xmin) into 100 bins, wc Computed the probability 
for each value of x to belong to the bin i (1 < i < 100). 
The resulting probability distribution P{xi/x) (normal- 
ized to unity) is drawn in Fig. versus the ratio of the 
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average susceptibility Xi of the bin i and the global aver- 
age susceptibility x- 




Fig. 5. Probability distribution of the susceptibility versus the 
ratio of the average value of the susceptibility of bin i and the 
global average susceptibility x for the bond-diluted Ising model 
(p = 0.7, 0.55, and 0.4) for L = 40, 64, and 96. The simulations 
are performed at Kc{L). The vertical dashed line indicates the 
average susceptibility Xi/x ~ 1- 



For the three studied dilutions, the shape of the curve 
does not become sharper with increasing size and the gen- 
eral shape remains remarkably independent of the dilu- 
tion, a phenomenon which strongly contrasts with the sit- 
uation encountered in the case of the q = A Potts model ■ 
Non self-averaging can be quantitatively checked by evalu- 

2 

ating the normalized squared width R^iL) = V-^(L)I x{L) , 
where is the variance of the susceptibility distribu- 

2 

tion: V^{L) = x^{L) — x(^) • The same quantity is also 

2 

evaluated for the magnetisation: Rm{L) — Vm{L)/m{L) . 
These ratios are shown versus the inverse lattice size for 
the three studied concentrations of the disordered Ising 
model in Fig. El The fact that Rm and R^ approach a con- 
stant when L increases as predicted by Aharony and Har- 
ris [71] is the signature of a non-self-averaging system, in 
agreement with the results of Wiseman and Domany |43l 
1^ for the site-diluted 3D Ising model. In fact, our numer- 
ical values of Rm are close to the estimate R^ = 0.055(2) 
for site dilution reported in Ref. when the system 
size goes to infinity, our data are compatible with a con- 
vergence towards Rm = 0.05(1), where the error estimate 
reflects the uncertainty in the remaining finite-size correc- 
tions which appear to be more pronounced than for site 
dilution. On the other hand, our limiting value for R^ 
seems at first sight to be completely off from the estimate 
R^i = 0.156(4) given in Ref. Upon closer inspection 
this discrepancy can be traced back to different definitions 
of the susceptibility. While Wiseman and Domany worked 
with the "high-temperature" expression x' = KV{ii^), 




0.00' ' ' ' ' ' ' ' ' ' ' ' 1 

0.02 0.04 0.06 0.08 0.1 0.12 

1/L 

Fig. 6. Normalized squared width of the magnetisation and 
susceptibility distributions versus the inverse of the lattice size 
for the three concentrations p — 0.4, 0.55, and 0.7 at the effec- 
tive critical coupling Kc{L). The straight lines are linear fits 
used as guides to the eye. 



we used the connected correlator x = KV{{^'^) — (|/i|)^), 
valid in both the low- and high-temperature phase. An- 
other difference is that they evaluated x' at the infinite- 
volume critical coupling K"^ while we followed the pseudo- 
transition points -ft'maxl.^) defined from the location of the 
X maxima (see below). In fact, when we evaluate R^i (i.e., 
without subtraction), we find for p = 0.7 and the largest 
sizes (where if max (i) ~ K"^) about 10 times larger values 
than shown in Fig. E| i?^' = 0.1635, 0.1766, and 0.1630 
for L — 64, 80, and 96. This is in accord with a remark 
in Ref. who noticed also for site dilution that R^ is 
smaller by a factor of 7—10, which would lead to a crass 
violation of the leading-order prediction in e = 4— d for the 
"ratio of ratios" Rm/Rx = 1/4 |7H. while working with y' 
they obtained a much closer ratio of Rm/Rx' = 0.35(2). 
In our p = 0.7 case we get, with Rm = 0.0407, 0.0423, 
and Rm — 0.0399 for L = 64, 80, and 96, estimates of 
Rm/Rx' = 0.2489, 0.2395, and 0.2448, in perhaps sur- 
prisingly good agreement with the renormalization group 
prediction of 1/4. Of course, given the large corrections 
seen otherwise, this agreement may be accidental, but it 
is also conceivable that the corrections do indeed partially 
cancel in the "ratio of ratios" . 

Still, the apparent non-universality of R^ and to a 
weaker extent also of Rm (the ratios seem to depend on 
the impurity concentration) remains puzzling. It may be 
attributed to crossover effects between the pure, disorder 
and percolation fixed points which may be more sensitive 
for bond dilution than for site dilution. Another possi- 
ble reason for quantitative discrepancies might be the fact 
that Rm and R^ were estimated in our study at the tem- 
perature corresponding to the maximum of the average 
susceptibility for the corresponding lattice size, and not 
at the critical point of the infinite system. This may lead 
to different scaling limits (similar to the different Binder 
parameter scaling limits for pure systems in the vicinity 
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Table 2. Values of the simulation points K = J/ksT (ex- 
tremely close to the susceptibility maxima) for the three mainly 
studied dilutions. 



L 


p = 0.70 


p = 0.55 


p = 0.40 


10 


0.6560 


0.8680 


1.3500 


12 


0.6556 


0.8680 


1.3210 


14 


0.6550 


0.8675 


1.3185 


16 


0.6546 


0.8655 


1.3170 


18 


0.6546 


0.8655 


1.3175 


22 


0.6542 


0.8655 


1.3160 


26 


0.6541 


0.8650 


1.3147 


30 


0.6538 


0.8650 


1.3144 


35 


0.6538 


0.8650 


1.3142 


40 


0.6538 


0.8655 


1.3141 


50 


0.6537 


0.8653 


1.3142 


64 


0.6535 


0.8655 


1.3144 


80 


0.6535 


0.8649 


1.3136 


96 


0.6535 


0.8649 


1.3136 



1500 




Fig. 7. Variation of the magnetic susceptibility and of the 
magnetisation (in the insert) versus the coupling strength K = 
J/ksT for the bond-diluted Ising model (p = 0.7) and the 
lattice sizes compiled in Table |21 going from 10, 12, 14, . ..to 
96. For each value of L, the curves are obtained by standard 
histogram reweighting of the simulation data at the value of 
K given in Tabled 

of Tc), albeit causing presumably only rather small devi- 
ations. 



4 Critical behaviour 

In order to study the critical behaviour of the bond-diluted 
Ising model, we have concentrated on the three particular 
concentrations p = 0.7, 0.55, and 0.4, for which the simu- 
lated lattice sizes go up to L = 96. For interested readers, 
the values of the simulation temperatures are reported in 
Table 



In Fig. Owe show results of typical runs for the mag- 
netisation (in the insert) and the susceptibility for p = 0.7 
at a particular temperature (or coupling) very close to the 
size-dependent critical point using the histogram reweight- 
ing technique. From the scaling of the location and height 
of the peaks with lattice size of this and similar quantities 
we then extracted the critical exponents of the system. 



4.1 Finite-Size Scaling Study 

4.1.1 Correlation length exponent 

As mentioned in the introduction, the vicinity of the ex- 
ponent ratios P/v and ^/v for the pure and disordered 
universality classes does not allow, by the use of stan- 
dard finite-size scaling (FSS) techniques, to discriminate 
between the two fixed points from the behaviour of the 
magnetisation and the susceptibility, respectively. Only 
the critical exponent which can be evaluated from the 
asymptotic FSS behaviour of the derivative of the mag- 
netisation versus temperature, 

d\nm/dK ^ adinniL^^" , (5) 

will be useful with this technique. But even for v, one is 
trying to resolve an expected shift from the pure model's 
value by less than 10%. Taking the data evaluated at 
Kraax and assuming the leading power-law behaviour Q , 
we have extracted by (linear) two-parameter least-squares 
fits over successively smaller ranges Lmin — Lmax the 
size-dependent effective exponent {l/v)eff- Using the first 
set of data for lattice sizes L = 4, 6, 8, . . . , 20 = imax, the 
resulting exponents are plotted in a broad range of concen- 
trations p against 1/L,„in in Fig. [S] (L,nin is the smallest 
lattice size used in the fits). We see that in the regime 
of low dilution {p close to 1), the system is clearly influ- 
enced by the pure fixed point. On the other hand, when 
the bond concentration is small, the vicinity of the per- 
colation fixed point induces a decrease of 1/v below its 
expected disordered value. Indeed, the percolation fixed 
point is characterized by l/i^ « 1.12 66 . 

As is demonstrated in Fig. IHl the same analysis for 
the three mainly studied dilutions with considerably larger 
lattice sizes going up to i = imax = 96 confirms the ob- 
servation of the great influence of both the pure and per- 
colation flxed points. The left-most data points in Fig. |51 
with Lmin = 40 follow from two-parameter least-squares 
fits of {d\nfh/dK)K^^^ for our five largest lattice sizes 
L = 40, 50, 64, 80, and 96. Since this leaves only three 
independent degrees of freedom, moving further into the 
asymptotic large-L regime would be meaningless with the 
data available to us. The variation of (l/i^)cff with L,„in 
clearly indicates deviations of the data from the asymp- 
totic ansatz rooted in confluent corrections and cross- 
over terms or both. As will be discussed in more detail be- 
low, we found it impossible, however, to resolve these cor- 
rections within a non-linear four-parameter ansatz. Still, 
some rough estimates of the critical exponent v can be 
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read off from Fig. IHl which are collected in Table 13 Ob- 
viously, it would be very difficult to extract more pre- 
cise values from this analysis. At a qualitative level, how- 
ever. Fig. El clearly indicates that the dilution for which 
the cross-over influence will be the least is around p ~ 
0.55 which suggests that the scaling corrections should be 
rather small for this specific dilution. 



4.1.2 Critical couplings 

The peak locations of the average susceptibility determine 
with a good accuracy the size-dependent (pseudo-) criti- 
cal couplings Kc{L), and from an infinite-size extrapola- 
tion, Kc{L) = + axL'^^'^ , the critical coupling in the 
thermodynamic limit can be estimated. Inserting the just 
determined estimate of ly, the linear fit shown in Fig. 11111 
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1 1 1 , 


'1 




- 
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- 








oo 









,1 
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Fig. 10. Size-dependent critical coupling Kc{L) versus L"^^" 
for the bond-diluted Ising model with p = 0.7, using our esti- 
mate oi V = 0.66(1). 



for p = 0.7 yields: 



= 0.6534(1). 



(6) 



The same procedure gives for p = 0.55 and p = 0.4 the 
critical couphngs = 0.8645(2) and = 1.3129(3), 
respectively. 



4.1.3 Susceptibility and magnetisation exponents 

As already mentioned, the average magnetisation to and 
susceptibility x are expected to scale at the size-dependent 
critical coupling Kc{L) with the lattice size as: 



1,-^ 1 



Xn 



(7) 



where a„i and are non- universal amplitudes. From least- 
squares fits to our extensive Monte Carlo data, we com- 
puted the effective size-dependent ratios of the critical ex- 
ponents {(3/v)cs and (7/i^)cff which are plotted in Fig. 1111 
against 1/Lmin for p = 0.7, 0.55, and 0.4, where imin is 
again the smallest lattice size used in the fits. 

Concerning the magnetisation, the effective ratios clear- 
ly converge, when the size increases, towards 0.515(5) for 
the three dilutions in agreement with the expected values 
of the pure and site-diluted models. For the susceptibil- 
ity, on the other hand, the behaviour of (7/1^)05 is more 
differentiated depending on the concentration p: the cases 
p — 0.7 and 0.55 are compatible with 1.965(10) if we take 
into account the error bars, very close to the pure and 
site-diluted model values, but the case p = 0.4 displays 
a (7/i^)off exponent a little bit larger. This discrepancy 
is probably due to the influence of the percolation fixed 
point whose ^/v ratio is close to 2.05, according to the 
same scenario as for the exponent 1/v. The critical expo- 
nents from the FSS study are summarized in Table |3| for 
the three studied dilutions, in good agreement with the 
hyperscaling relation d = 2(3 / v -\- ^ /v. 

Following the work of Ballesteros et al. 021 which em- 
phasized the great influence of corrections-to-scaling to get 
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Table 3. Critical exponents deduced from the FSS study of 
the three dilutions. 



p 


0.7 


0.55 


0.4 




1.52(2) 


1.46(2) 


1.42(2) 


V 


0.660(10) 


0.685(10) 


0.705(10) 


13/u 


0.515(5) 


0.513(5) 


0.510(5) 




1.965(10) 


1.977(10) 


2.000(10) 


213 /v + -f/u 


2.995(20) 


3.003(20) 


3.020(20) 



an accurate determination of the critical exponents in the 
site-diluted Ising model, we tried to estimate them in the 
bond-diluted version as well. If we consider an irrelevant 
scaling field g with scaling dimension yg — —to < 0, the 
scaling expression of, e.g., the susceptibility is: 



x(£-^^,ff) = i^/Vx(i^^i"".9)■ 



(8) 



At t = 0, around the fixed point value 5 = 0, this leads to 
the standard expression a^L^/^ll+h^L-^ ^OiL-"^"^)]. We 
hence tried to fit the susceptibility data using the ansatz: 



(9) 



We systematically varied the exponent ratio ^/v and the 
corrections-to-scaling exponent uj and determined the min- 
imum of x^/d.o.f by performing linear fits in the range 
18 < L < 96 to fix the amplitudes and h^. The result- 
ing x^-landscape for p = 0.7 is shown in Fig. 1121 where 
the base plane was restricted to the range 1.8 <^/v< 2.2 
and < w < 4. The absolute minimum was found at 
"f/v = 1.952, ijj — 2.24, and a second, less pronounced 
minimum at 7/t/ = 2.112, oj = 0.01. The results for the 
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11. Effective ratios {'y/i')aB and [j3/i')e« as a function of 
min for p — 0.7, 0.55, and 0.4. The error bars correspond 
to the standard deviations of the power-law fits. The arrows 
indicate roughly the expected values for both the pure and 
disordered model. 



Fig 

1/L 



magnetisation respectively its logaritmic derivative and 
the other dilutions p = 0.55 and 0.4 look qualitatively 
similar. 

The figure may again be interpreted in favor of a com- 
petition between the two fixed points: the first one is char- 
acterized by a valley in the w-direction, with an absolute 
minimum for a large value of uj, i.e., small corrections-to- 
scaling since these corrections scale as L~". Then, this 
fixed point should be the disordered one and the second 
the percolation fixed point for which "f/v ~ 2.05 and 
P/v ~ 0.48 iS^. But the accurate determination of the 
corrections-to-scaling exponent to is very difficult because 
it strongly depends on the quantity considered (xmax, 
"^ifmax oi' ('^l'^''^/'^^)^max) ^-ud ou thc sizcs iucludcd in 
the fits. Indeed, the valley in the uj direction does not dis- 
play any deep minimum and the absolute minimum ob- 
served does not seem to be relevant concerning the value 
of to. We are thus not able to extract a numerical value of 
the corrections-to-scaling exponent, in contrast to previ- 
ous claims that uj w 0.4 |I^48U49| . 



4.2 Temperature Scaling 

The temperature behaviour of the susceptibility and mag- 
netisation should allow to avoid the previous difficulties 
generated by the vicinity of the ratios of the critical ex- 
ponents. The critical exponents are indeed the following: 
pure Ising model: /3 = 0.3258(14), 7 = 1.2396(13) |S2|, 
site-diluted model: f3 = 0.3546(28), 7 = 1.342(10) gl]. 

As a function of the reduced temperature t — (Kc — K) 
(i < in the low-temperature (LT) phase and t > in the 
high-temperature (HT) phase) and the system size L, the 
magnetisation and the susceptibility are expected to scale 




Fig. 12. Plot of the x^-landscape deduced from linear fits of 
Xmax(i) = a^L'^^'^il + h^L-'^) in the range 18 < L < 96 
for the bond-diluted Ising model with concentration p — 0.7, 
exhibiting an absolute minimum at y/v = 1.952, u) = 2.24 
and a secondary minimum at j/i/ — 2.112, uj = 0.01. The 
exponents are fixed parameters and the amplitudes are free. 
The cutoff at ~ 300 has been introduced in order to improve 
the clarity of the figure. 
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Fig. 13. Variation of the temperature dependent effective crit- 
ical exponent /3eff(|t|) ~ dlnm/dln \t\ (in the low-temperature 
phase) as a function of the reduced temperature \t\ (top) and 
L^^'^\t\ (bottom) for the bond-diluted Ising model with p = 0.7 
and several lattice sizes L. The data points are directly ob- 
tained from Monte Carlo data (no histogram reweighting). The 
magnetisation vs. the coupling strength K = J/ksT in the or- 
dered phase K > is shown in the upper part. 

as: 

mit,L)^\tff±iL'/''\t\), (10) 

x{t,L)^\trg±iL'/''\t\), (11) 

where f± and g± are scaling functions of the variable 
X = L^/^lt] and the subscript ± stands for the HT/LT 
phases. Then we can define temperature dependent effec- 
tive critical exponents 

P,ei\t\)=dlnm/d\n\t\, (12) 
7cff(|i|) = -dlnx/rfln|t|, (13) 

which should converge towards the asymptotic critical ex- 
ponents /3 and 7 when L 00 and \t\ — > 0. The results 
for p = 0.7 are shown in Figs. 1131 and 1141 

The effective exponents /3eff(|i|) evolve as a function 
of |t| between a maximum value and in the two |t| di- 
rections: when |t| is large {K far from Kc), the system 
is outside the critical region and the description with the 
critical exponents is no more valid whereas when \t\ —^ 0, 
the correlation length becomes of the same order as the 
linear size of the system and the finite-size effects be- 
come very strong. Nevertheless, just before the curves 
bend down, for the largest lattice sizes we can read off that 
i^eff (1^1) ~ 0.34 — 0.36. In the case of the susceptibility, for 
the greatest sizes, 7cff(|i|) is stable around 1.34 when |t| 
is not too small, i.e., when the finite-size effects are not 
too strong. The plot of 7cff(|t|) vs. the rescaled variable 



Fig. 14. Variation of the temperature dependent effective criti- 
cal exponent 7cff (|t|) ~ — dlnx/rfln \t\ (in the low-temperature 
phase) as a function of the reduced temperature \t\ (top) and 
L^^^\t\ (bottom) for the bond-diluted Ising model with p = 0.7 
and several lattice sizes L. The horizontal dashed lines indicate 
the pure and site-diluted values of 7. The susceptibity vs. the 
coupling strength K = J/ksT in the ordered phase is shown 
in the upper part. 



L^/'^\t\ shows that the critical power-law behaviour holds 
in different temperature ranges for the different sizes stud- 
ied. As expected, the size effects are more sensitive when 
the lattice size is small and the critical behaviour is better 
described when the size increases. 

A more quantitative analysis is possible for the con- 
centrations p — 0.55 and 0.7 for which the most accurate 
temperature simulations have been done. /^From the tem- 
perature behaviour of the susceptibility, we have directly 
extracted the power-law exponent 7 from error weighted 
least-squares fits by choosing the temperature range that 
gives the smallest x^/d-O.f for several system sizes. The 
results are given in Tabled We see that all estimates are 
consistent with 7 « 1.34 — 1.36, which is clearly different 
from the pure model's exponent of 7 « 1.24. In particular 
we do not observe in this analysis any pronunced differ- 
ences between the two dilutions. For p — 0.4, our main 
dilution closest to the percolation threshold, the number 
of accurate temperature points is unfortunately not large 
enough to allow for reliable fits. 

From the previous expressions of the magnetisation 
and susceptibility as a function of the reduced temper- 
ature and size, following a procedure proposed by Binder 
and Landau |72j a long time ago, it is instructive to plot 
the scaling functions f±{x) and g±{x). For finite size and 
\t\ ^ 0, the scaling functions may be Taylor expanded in 
powers of the inverse scaling variable x'^ — {L^^'^\t\)^^ , 
for example in the case of the susceptibility x±(^5^) — 
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Fig. 15. Log-log plot of the scaling functions mL^^" (top) 
and xi"^''" (bottom) against L^^''\t\ for p = 0.4, 0.55, and 0.7 
in the low-temperature ordered phase {K > K^). The lines 
show the power-law behaviours with the exponents (3 ~ 0.355 
and 7 ~ 1.34 which characterize the disorder fixed point. The 
data for the smallest values of L^^'^\t\ which do not lye on the 
master curve show the size effects described previously when 
the correlation length ^ is limited by the linear size L. 

Table 4. Critical exponent 7 deduced from log-log fits of the 
susceptibility vs. the reduced temperature \t\ in the ordered 
phase for the concentrations p = 0.55 and 0.7. 



p = 0.55 L 


40 




22 




# points 


7 




12 




xVd.o.f 


0.07 




1.6 




7 


1.36(3) 




1.36(2) 




p = 0.7 L 


40 


35 


22 


18 


# points 


10 


8 


9 


8 


xVd.o.f 


8.37 


0.01 


0.08 


0.81 


7 


1.32(1) 


1.34(5) 


1.40(4) 


1.34(5) 



\t\ '^[g±{'X')+x g'^{oo) + 0{x ^)], where the amplitude 
g±{oo) is usually denoted by r±. Multiplying by L^''^" 
leads to 



X± 



(14) 



where ^±(2^) = x~'^ g±(x). The case of the magnetisation 
is slightly different, since the magnetisation is asymptot- 
ically vanishing in the high-temperature phase, and thus 
rh±{t,L) = |t|^[/±(oo) +a;-i/4(oo) -F0(x-2)] 







B'+x-^ 
B'_x-^ 



0(x-2) 
0(x-2^ 



(15) 

which should give universal curves for the different sizes 
and temperatures. 



Fig. 16. Log-log plot of the scaling functions mL"/'^ (top) 
and xL~''^'^ (bottom) against L^^'^\t\ in the disordered phase 
{K < K^) for the three dilutions. The lines show the power- 
law behaviours with the exponents /3 — 1 ~ —0.645 and —7 ~ 
— 1.34 which characterize the disorder fixed point. 

The curves in the ordered phase shown in Fig. [T^ 
are obviously universal master curves whose slopes, in a 
log-log plot, give the critical exponents (3 ~ 0.355 and 
7 ~ 1.34. Indeed, when \t\ — > but with L still larger 
than the correlation length ^, one should recover the crit- 
ical behaviour given by /-(a;) ^ x^ and g~{x) ^ x~^ . 
The same procedure applied to the two other bond con- 
centrations p = 0.55 and 0.4 gives analogous results and 
confirms the universal values of the disordered critical ex- 
ponents /3 and 7, independent of the concentration as is 
illustrated in Fig. El 

In the disordered phase K < K^, t > the previ- 
ous scaling assumptions for fh and x l^ad to universal 
functions whose slopes on a logarithmic scale are equal to 
/3 - 1 ~ -0.645 and -7 ~ -1.34, see Fig-UHl 

Some combinations of critical amplitudes are also uni- 
versal and thus characterize the critical point. In fact, crit- 
ical amplitude ratios are often more sensitive to the uni- 
versality class than the critical exponents themselves |73] . 
Among such ratios, F^/F^ is directly accessible through 
our results. An example is illustrated in Fig. 1171 where the 
scaling function of Eq. pi|l is plotted against the scaling 
variable x. The values of the amplitudes which describe 
the approach to criticality from above and from below, 
F^ and respectively, are shown in shaded horizontal 
stripes. The ratios 

= 1.62 ±0.10 (p = 0.7), (16) 
= 1.50 ±0.10 (p = 0.55), (17) 
F+/F^ = 1.48 ± 0.20 {p = 0.4) (18) 

follow. The values obtained for the three dilutions are con- 
sistent within error bars but they unfortunately appear 
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Fig. 17. Log-log plot of the scaling functions xl^T against 
L^^'^\t\ for the three dilutions. The shaded horizontal stripes 
indicate the critical amplitudes /+ and r-{< / + ), respectively. 

to be in contradiction with a field-theoretic approach of 
Bervillier and Shpot 60 , who predicted a ratio equal to 
r^/F- = 3.05(32). This is nevertheless not a crucially 
conflicting result, since susceptibility amplitudes of course 
depend on the definition used to compute the susceptibil- 
ities. At any rate, our estimates (|16|I ~ H18() for the disor- 
dered Ising model are clearly different from those for the 
pure model where / + / F- obtained with field-theoretic 
and high-temperature series expansion techniques varies 
between 4.70 and 4.95 [Tl], and a recent high-precision 
Monte Carlo study |7S| concluded that F+/F^ = 4.75(3). 



5 Conclusions 

In this paper, we have reported an intensive Monte Carlo 
study of the physical properties of the 3D disordered bond- 
diluted Ising model. As the critical exponents of the pure 
and disordered models are very close, the numerical pro- 
cedure has to be carefully chosen, especially the thermal 
averaging and the average over the disorder realisations. 
We have determined numerically the phase diagram of the 
disordered model, in very good agreement with a single- 
bond effective-medium approximation. 

The non-self-averaging at criticality, as discussed by 
Aharony and Harris, is confirmed here. At a quantita- 
tive level, the discrepancy between the numerical values of 
the normalized squared width of the susceptibility R^{L) 
in this work and that in Ref. 0^ definitely comes from 
the difference in the definition of x- Although in our case 
Rx{L) and to a lesser extent also Rm{L) do not seem to 
approach a unique constant independent of the disorder 
concentration when the system size L increases, we cannot 
exclude a unique limiting value when cross-over effects are 
properly taken into account. This possibility is corrobo- 



Table 5. Critical exponents and critical amplitude ratio of the 
susceptibility for the different fixed points. 



fixed point /3 7 



fixed point P/i^ 'y/p F+fF^ 



pure 


0.517(3) 


1.966(6) 


4.75(3) 


percolation 


0.45 


2.05 




site-diluted 


0.519(3) 


1.963(5) 




bond-diluted 


0.515(5) 


1.97(2) 


1.50(20) 



rated by our finite-size scaling study of critical exponents. 
A comparison of the underlying disorder distributions for 
the bond-diluted Ising and q = 4 Potts models shows that 
the weight of rare events plays a much more crucial role 
for the critical properties of the latter model. 

The main emphasis was on investigations of the crit- 
ical behaviour of the disordered model based on a finite- 
size scaling study and on a complementary study of the 
temperature scaling behaviour. In the first case, we have 
shown the great infiuence of the cross-over phenomena 
between the pure, disorder and percolation fixed points 
which lead to the measurement of effective critical expo- 
nents dependent on the concentration of magnetic bonds. 
The temperature study, on the other hand, turned out 
to be much better behaved and allowed us to confirm the 
unique values of the critical exponents, independent of the 
concentration and in agreement with the site-diluted val- 
ues, thus providing strong numerical evidence for the theo- 
retically predicted universal behaviour of the two disorder 
distributions. These results are summarized in Table [S] A 
scaling study of the average of the maxima of the suscepti- 
bilities Xj for each sample instead of the maximum of the 
average susceptibility would have been an interesting al- 
ternative. However, it is technically more complicated be- 
cause from sample to sample the position of the maximum 
of the susceptibility varies and may be far from the simu- 
lation temperature making histogram reweighting unreli- 
able. There should even appear samples with large bond- 
density fluctuations where unconnected or weakly con- 
nected clusters of bonds coexist. Such special configura- 
tions may display a double-peaked (or more complicated) 
susceptibility curve resulting from the independence of 
these uncorrelated clusters. 

Finally, we should remark that even though this work 
required a really huge amount of effort and CPU time, 
for some reasons that we are unable to identify, the ac- 
curacy of our results is not as good as in the work of 
Ballesteros et al. [23 for site-dilution. For example we 
were not able to conclude for a reliable estimate of the 
corrections-to-scaling exponent w as they did. Neverthe- 
less it is worth noticing that the temperature scaling is 
quite satisfying and led to the measurement of the suscep- 



pure 0.6304(13) 0.3258(14) 1.2396(13) 

percolation 0.89 0.40 1.82 

site-diluted 0.6837(53) 0.3546(28) 1.342(10) 

bond-diluted 0.68(2) 0.35(1) 1.34(1) 
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tibility amplitude ratio I+Z-T-, a universal combination 
which was still unknown numerically, and which comple- 
ments an estimate of other universal combinations of am- 
plitudes reported recently in Refs. 03211 ■ As anticipated 
from a previous analysis for disordered systems, our esti- 
mates for the amplitude ratio /+ / are clearly different 
from the pure model and thus yield the cleanest signal for 
a unique disorder fixed point. 
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